Bienvenid@s a la uuuuultima tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la ultima parte del curso, los cuales se enfocan principalmente en aplicar inferencia bayesiana para generar regresiones lineales y estudiar métodos de obtención de la posterior mas poderosos, como es MCMC. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Videos de las clases:
Documentación:
Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.
Cross-Validation
La cross-validation divide los datos de entrenamiento en
\(k\) subconjuntos, entrenando el
modelo \(k\) veces, siempre dejando uno
de los subconjuntos fuera. Esto logra que el modelo no se ajuste a unos
datos en particular, sino que aprenda de mejor manera de todos los
datos.
Criterios de Información
Los criterios de información permiten, de cierta forma,
medir cuanto overfitting o underfitting tienen distintos modelos
entrenados con los mismos datos. Esto permite comparar distintos modelos
para ver cual es el que se comporta mejor. Aún así, esto no resuelve el
problema de overfitting, solo lo mide, ya que estas medidas mejoran en
general mientras más complejo (y por lo tanto, más ajustado) es el
modelo a los datos.
Regularización
La regularización permite mantener los parámetros del
modelo controlados al agregar una penalización en la función objetivo a
qué tan grande sean estos parámetros. Esto es controlado por una
constante \(\lambda \in [0,1]\) que
controla la cantidad de overfitting del modelo. En general, penalizar
tanto los parámetros genera un riesgo de underfitting, ya que no permite
que los parámetros crezcan lo suficiente, por eso se suele ajustar
usando datos de validación.
Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.
Respuesta
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
# Manipulación de varios plots en una imagen.
library(gridExtra)
# Análisis bayesiano
library("rethinking")
Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:
install.packages("rethinking")
En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
install.packages(c("mvtnorm","loo","coda","remotes","shape"), repos="https://cloud.r-project.org/",dependencies=TRUE)
remotes::install_github("rmcelreath/rethinking")
El objetivo de esta pregunta es introducirlo en la aplicación de una
regresión bayesiana. Con esto, buscaremos entender como calcular una
regresión bayesiana en base al “motor” aproximación de Laplace,
revisando como se comporta la credibilidad de sus predicciones y como la
regresión lineal puede llegar a mutar a aplicar una transformación en el
vector \(x\). Para responder esta
pregunta centre su desarrollo solo en las clases y las funciones que nos
brinda la librería rethinking.
Unos expertos en alometría deciden realizar un estudio de las alturas
de unos niños en un colegio, buscando generar un regresor lineal
bayesiano capaz de predecir la altura en base al peso de los alumnos.
Para realizar este trabajo recopilan los datos
table_height.csv, quien posee observaciones fisiológicas de
192 alumnos.
Parte I
En conocimiento los datos recopilados por los expertos, le solicitan realizar la siguiente serie de tareas:
Parte II
En base a los resultados obtenidos, el experto que trabaja con usted le señala que las alturas se suelen modelas con pesos logarítmicos, por lo que le sugiere añadir un logaritmo natural en el vector \(x\) que compone su modelo lineal. Realice nuevamente la regresión utilizando un intervalo del \(95\%\) sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos, señalando si el modelo logra ajustar mejor los valores.
Respuesta:
Respuesta Aquí
Bonus
Esta parte es opcional y contará como un puntaje adicional si la realizan, pero si no lo hacen no habrá penalización.
Compare los resultados obtenidos con el método de Laplace con MCMC y comente los resultados. Para esto pueden utilizar ulam del paquete rethinking.
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) #necesitarán instalar estas librerías
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# si no les funciona instalar cmdstan con lo anterior pueden probar con:
check_cmdstan_toolchain(fix = TRUE)
install_cmdstan()
Respuesta Aquí
El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.
En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:
\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \]
De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).
Respuesta:
Primero, definimos la función gamma_func para calcular
los valores de la función de densidad de la distribución gamma, según la
definición dada anteriormente \[
f(x\mid \alpha,\beta) =
\begin{cases}
x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\
0 ~&\text{si } x \leq 0
\end{cases}
\]
gamma_func <- function(x, shape, rate) {
if (x > 0)
return(x^(shape - 1) * exp(-rate * x))
else
return(0)
}
Luego, definimos la función metropolis la cual nos va a
permitir calcular lo pedido
metropolis <- function(theta_0, n_iter, alpha, beta) {
theta <- vector(length = n_iter + 1)
current <- theta_0
for (k in 1:n_iter + 1) {
theta[k] <- current
# Obtenemos el valor propuesto usando la distribución normal
proposal <- rnorm(1, mean = theta[k], sd = 1)
# Valor de la densidad gamma para este valor propuesto
prob_prop <- gamma_func(proposal, shape = alpha, rate = beta)
# Probabilidad para el valor actual
prob_current <- gamma_func(current, shape = alpha, rate = beta)
# Proporción entre ambas probabilidades
ratio <- prob_prop / prob_current
# Probabilidad de aceptar el nuevo valor
prob <- min(1, ratio)
# Decisión sobre si aceptar el nuevo valor o no
decision <- rbinom(1, 1, prob)
# Actualización del valor
current <- ifelse(decision == 1, proposal, current)
}
# Devolvemos theta_1 a theta_(n+1), dando n valores
return(theta[-1])
}
Ahora, graficamos el histograma para la cantidad de valores \(10^1\), \(10^2\), \(10^3\), \(10^4\), \(10^5\) y \(10^6\), considerando \(\theta_0 = 1\)
hist(metropolis(1, 1e1, 5, 0.2))
hist(metropolis(1, 1e2, 5, 0.2))
hist(metropolis(1, 1e3, 5, 0.2))
hist(metropolis(1, 1e4, 5, 0.2))
hist(metropolis(1, 1e5, 5, 0.2))
hist(metropolis(1, 1e6, 5, 0.2))
Notemos que a partir de las 1000 iteraciones, el gráfico se comienza a parecer a la distribución gamma. Mientras más iteraciones, mejor se ajusta el modelo.
Ahora cambiemos la condición inicial. Usaremos los valores \(\theta \in \{1, 5, 10, 50, 100\}\) con \(n_{iter} = 5 * 10^5\)
hist(metropolis(1, 5e5, 5, 0.2))
hist(metropolis(5, 5e5, 5, 0.2))
hist(metropolis(10, 5e5, 5, 0.2))
hist(metropolis(50, 5e5, 5, 0.2))
hist(metropolis(100, 5e5, 5, 0.2))
Notemos que en general la forma de los gráficos no cambia dependiendo
de la condición inicial. El algoritmo MCMC converge
igualmente a la distribución buscada. Lo que podría pasar es que, si la
condición inicial es “muy mala”, el algoritmo necesite más iteraciones
para converger.
Ahora, utilizando \(\theta_0 = 1\) y
\(n_{iter} = 10^6\), comparemos el
gráfico con la función gamma nativa de R,
rgamma.
Para esto, graficaremos la función de densidad gamma usando la
función rgamma. Para esto, generaremos \(10^6\) elementos aleatorios considerando
\(\alpha = 5\) y \(\beta = \frac{1}{5} = 0.2\)
# Distribución no normalizada
gamma_vect <- rgamma(1e6, shape = 5, rate = 0.2)
hist(gamma_vect, xlim = range(c(0, 120)))
Luego, grafiquemos usando la función metropolis
definida
# Distribución no normalizada
gamma_vect <- metropolis(1, 1e6, 5, 0.2)
hist(gamma_vect, xlim = range(c(0, 120)))
Notemos que, efectivamente, ambas gráficas son bastante similares en
forma, lo que nos permite afirmar que el algoritmo
metropolis se ajusta bien a la distribución
gamma y por lo tanto permite estimarla correctamente.
A work by CC6104